Playing Around with Fitness Functions

Imports

Global Variables

popSize <- 500
nChr <- 10
nSegSites <- 100
nGens <- 50

Functions

oneTraitFitFunc <- function(x) {
  return (-(x)^2) + rnorm(1,sd=2)
}

twoTraitFitFunc <- function(x) {
  res = -(x[,1]^2) - x[,2]^2
  return (res)
}

calculateFitnessTwoTrait <- function(x,y) {
  return ((-(x)^2)+(-(y)^2))
}

hetLocus <- function(locus) {
  return (length(unique(locus)) > 1)
}

Plotting Functions

theme <- theme(plot.background = ggplot2::element_blank(),
               panel.background = ggplot2::element_blank(),
               axis.line = ggplot2::element_line(linewidth = 0.2),
               plot.title = ggplot2::element_text(hjust = 0.5,
                                                  face = 'bold',
                                                  size = 12),
               axis.text = ggplot2::element_text(size  = 12,
                                                 color = 'black'),
               axis.title = ggplot2::element_text(size  = 12,
                                                  face = 'bold',
                                                  color = 'black'))

# Expects a dataframe with 2 columns:
# traitValA, traitValB (for type "CONTOUR"),
# with a 3rd column (fitness) (for type "SURFACE")
overlayWalkOnLandscape <-function(df,
                                  type="CONTOUR",
                                  fitCalc,
                                  traitAMin=-10,
                                  traitAMax=10,
                                  traitBMin=-10,
                                  traitBMax=10) {

  fitness_x = seq(traitAMin,traitAMax, by=0.25)
  fitness_y = seq(traitBMin,traitBMax, by=0.25)
  fitness_z = outer(fitness_x,fitness_y,fitCalc)
  fig <- plot_ly()

  if (type == "CONTOUR") {
    plot_ly() %>%
      layout(xaxis = list(title = "Trait A", constrain = "domain"),
                          yaxis = list(title = "Trait B", scaleanchor="x")) %>%
      add_trace(
        fig,
        x=fitness_x,
        y=fitness_y,
        z=fitness_z,
        type='contour',
        colorscale = "RdBu",
        colorbar=list(title = "w"),
        line = list(color = 'black', width = 1),
        opacity=1) %>%
      add_trace(
        fig,
        df,
        name = popSize,
        x = df$traitValA,
        y = df$traitValB,
        type='scatter',
        mode = 'lines',
        line = list(color = 'black', width = 4, dash = 'solid'),
        opacity = 1)
  } else if (type == "SURFACE"){
    plot_ly() %>%
      layout(scene = list(xaxis = list(title = "Trait A"),
                          yaxis = list(title = "Trait B"),
                          zaxis = list(title = "w"),
                          aspectmode='cube')) %>%
      add_trace(
        fig,
        df,
        name = popSize,
        x = df$traitValA,
        y = df$traitValB,
        z = df$fitness,
        type = 'scatter3d',
        mode = 'lines',
        opacity = 1,
        line = list(width = 10)
      ) %>%
      add_trace(
        fig,
        x=fitness_x,
        y=fitness_y,
        z=fitness_z,
        type='surface',
        colorbar=list(title = "w"),
        colors = "PuBuGn",
        opacity=0.9)
      
  } else {
    print("Type Not Supported")
  }
}

# Plot a population on a a 3d fitness chart
# only works for populations with 2 traits
# TODO: figure out how to add colorbar label
plot3dPopulationFitness <- function(pop, fitCalc) {
  popSize <- nInd(pop)
  df <- data.frame(traitA=numeric(popSize),
                   traitB=numeric(popSize),
                   fitness=numeric(popSize))
  for (i in 1:popSize) {
    df$traitA[i] <- pheno(pop)[i,1]
    df$traitB[i] <- pheno(pop)[i,2]
    df$fitness[i] <- fitCalc(df$traitA[i], df$traitB[i])
  }
  fig <- plot_ly()
  
  plot_ly() %>% 
    layout(scene = list(xaxis = list(title = "Trait 1"),
                        yaxis = list(title = "Trait 2"),
                        zaxis = list(title = "w"),
                        aspectmode='cube')) %>%
    add_trace(
      fig,
      df,
      name = popSize,
      x = df$traitA,
      y = df$traitB,
      z = df$fitness,
      type = 'scatter3d',
      mode = 'markers',
      color=df$fitness)
}

# Plots the trait architecture on a genetic basis
# Works for 1 or 2 trait populations
plotTraitArchitecture <- function(pop, fitFunc) {
  fit_pre <- mean(fitFunc(gv(pop)))
  geno <- pullSegSiteGeno(pop)
  cols <- colnames(geno)
  nLoci <- length(cols)
  popSize <- nrow(geno)
  eff_sizes <- data.frame(id = character(nLoci),
                          eff_size=numeric(nLoci))
  for (l in 1:nLoci) {
    id <- cols[l]
    eff_sizes$id[l] <- id
    locus = geno[,l]
    if (!hetLocus(locus)) {
      eff_sizes$eff_size[l] <- NA
      next
    }
    strs = unlist(strsplit(id, "_"))
    chr = strtoi(strs[1])
    site = strtoi(strs[2])
    pop1 <- editGenome(pop, ind=c(1:popSize), chr=chr, segSites=site, allele=0, simParam=SP)
    pop2 <- editGenome(pop, ind=c(1:popSize), chr=chr, segSites=site, allele=1, simParam=SP)
    effect_size_1 <- mean(fitFunc(gv(pop1))) - fit_pre
    effect_size_2 <- mean(fitFunc(gv(pop2))) - fit_pre
    if (effect_size_1 > 0) {
      eff_sizes$eff_size[l] <- effect_size_1
    } else if (effect_size_2 > 1) {
      eff_sizes$eff_size[l] <- effect_size_2
    } else {
      eff_sizes$eff_size[l] <- 0
    }
  }
  eff_sizes <- eff_sizes[apply(eff_sizes!=0, 1, all),]
  eff_sizes <- na.omit(eff_sizes)
  ggplot(data=eff_sizes, aes(x=reorder(id, -eff_size), y=eff_size)) +
    geom_bar(stat="identity") +
    labs(x = "Variant Id", y = "Effect Size") +
    theme +
    theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=1, size = 6, margin = margin(b = 10)),
          axis.text.y = element_text(margin = margin(l=10, r=10)))
}

# Show the fitness of a population
plot_fitness <- function(df) {
  ggplot(data=df, aes(x=gen, y=fitness)) +
    geom_line() +
    geom_point() +
    labs(x = "Generation", y = "Fitness")
}

# Plot Genetic Values for Two Traits
plot_hist <- function(pop) {
  gv_a = gv(pop)[,1]
  gv_b = gv(pop)[,2]
  idx = c(1:length(gv_a))
  df <- data.frame(gv_a, gv_b)
  df <- melt(as.data.table(df))
  ggplot(df, aes(x=value, color=variable)) + geom_histogram(binwidth=1, position='identity')
}

Play around with different fitness functions

x = seq(-100,100, by=1)
y = oneTraitFitFunc(x)

df <- data.frame(x=x, y=y)
ggplot(data=df, aes(x=x, y=y)) +
    geom_line() +
    geom_point() +
    labs(x = "Trait Value", y = "Fitness")

fitness_x = seq(-10,10, by=0.25)
fitness_y = seq(-10,10, by=0.25)
fitness_z = outer(fitness_x,fitness_y,calculateFitnessTwoTrait)


plot_ly(x=fitness_x,
        y=fitness_y,
        z=fitness_z,
        type='surface',
        colors = "PuBuGn",
        opacity=1) %>%
  layout(scene = list(xaxis = list(title = "Trait A"),
                      yaxis = list(title = "Trait B"),
                      zaxis = list(title = "w"),
                      aspectmode='cube'))

# Valid arguments:
# colors = "PuBuGn"
# colors = colorRampPalette(c("blue", "orange"))(15)
# colors = magma(50, alpha = 1, begin = 0, end = 1, direction = 1) (viridis, plasma, magma, inferno)

Show that all individuals converge to a fitness maximum

start_trait_val <- c(-50,-25,-10,10,25,50)
numSims <- length(start_trait_val)
popSize <- 500
nSegSites <- 50

df <- data.frame(gen=rep(1:nGens, times=numSims),
                 fitness=numeric(nGens*numSims),
                 initialVal = numeric(nGens*numSims),
                 traitVal=numeric(nGens*numSims))

for (s in 1:numSims) {
  initVal <- start_trait_val[s]
  founders = quickHaplo(
     nInd=popSize,
     nChr=nChr,
     segSites=nSegSites
  )
  SP <- SimParam$new(founders)
  SP$addTraitA(
    10,
    mean=initVal,
    var=abs(initVal/2)
  )
  SP$setVarE(H2=0.5)
  pop <- newPop(founders, simParam=SP)
  
  for(gen in 1:nGens) {
    idx = (s-1)*nGens + gen
    df$fitness[idx] <- oneTraitFitFunc(pheno(pop))
    df$traitVal[idx] <- meanP(pop)
    df$initialVal[idx] <- initVal
    pop <- selectCross(pop=pop, trait=oneTraitFitFunc, nInd=popSize/2, nCrosses=popSize)
  }
}

Warning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement lengthWarning: number of items to replace is not a multiple of replacement length

df$initialVal = as.character(df$initialVal)
ggplot(data=df, aes(x=gen, y=traitVal)) +
  geom_line(aes(color = initialVal)) +
  labs(x = "Generation", y = "Trait Value")

Plot Trait Value against Fitness

nGens <- 50
nInd <- 500
nSegSites <- 1000


founders = quickHaplo(
   nInd=nInd,
   nChr=nChr,
   segSites=nSegSites
)
SP <- SimParam$new(founders)
SP$addTraitA(
  10,
  mean=50,
  var=abs(initVal/2)
)
SP$setVarE(H2=0.5)
pop <- newPop(founders, simParam=SP)

fit_df <- data.frame(gen=1:nGens,
                 fitness=numeric(nGens),
                 traitVal=numeric(nGens))
for(gen in 1:nGens) {
  fit_df$fitness[gen] <- fitnessFunc(meanP(pop))
  fit_df$traitVal[gen] <- meanP(pop)
  pop <- selectCross(pop=pop, trait=oneTraitFitFunc, nInd=popSize/2, nCrosses=popSize)
}

par(mar = c(5, 5, 3, 5) + 0.3)
plot(fit_df$gen, fit_df$traitVal, type="l", lwd = "3", col=2, xlab = "Generation", ylab = "Trait Value")   
par(new = TRUE) 
plot(fit_df$gen, fit_df$fitness, type="l", lwd = "3", col = 3, axes = FALSE, xlab = "", ylab = "") 
axis(side = 4, at = pretty(range(fit_df$fitness)))
mtext("Fitness", side = 4, line = 3)
par(xpd=TRUE)
legend("right",
  c("Trait Value", "Fitness"),
       lty = 1,
       lwd = 3,
       col = 2:3)

Create 1 Trait Population

popSize = 200
nSegSites = 50
nLoci = nSegSites * nChr

founders = runMacs(
   nInd=popSize,
   nChr=nChr,
   segSites=nSegSites
)
SP <- SimParam$new(founders)
SP$addTraitA(
  mean=20,
  var=5,
  nQtlPerChr = 10 # changing this changes effect size? 100 QTL
)

SP$setVarE(H2=0.5)
foundingPop <- newPop(founders, simParam = SP)
plotTraitArchitecture(foundingPop, oneTraitFitFunc)

Create 2 Trait Population

set.seed(123)
popSize = 1000
nSegSites = 100

founders = runMacs(
   nInd=popSize,
   nChr=nChr,
   segSites=nSegSites
)
SP <- SimParam$new(founders)
SP$addTraitA(mean=c(runif(1,-10,10),runif(1,-10,10)), var=c(1,1), nQtlPerChr=10)

SP$setVarE(H2=c(0.5,0.5))
foundingPop <- newPop(founders, simParam = SP)

#plot3dPopulationFitness(foundingPop, calculateFitnessTwoTrait)
#plotTraitArchitecture(foundingPop, twoTraitFitFunc)

Out of alleles left in the population, how many of them would increase fitness? Next: how does this change with size of allele? Look at editGenomeTopQtl() Next: how does size of ‘fitness delta’ change over generations?

nGens = 20
pop <- foundingPop
selIntensity <- 0.2

cols = colnames(pullSegSiteGeno(pop))
df <- data.frame(gen=rep(1:nGens),
                 fitness=numeric(nGens),
                 traitVal=numeric(nGens),
                 percFitInc=numeric(nGens),
                 segAlleles=numeric(nGens))

for (gen in 1:nGens) {
  print(gen)
  df$fitness[gen] <- fitnessFunc(meanG(pop))
  df$traitVal[gen] <- meanG(pop)
  num_het_loci = 0
  num_inc = 0
  geno <- pullSegSiteGeno(pop)
  for (l in 1:nLoci) {
    strs = unlist(strsplit(cols[l], "_"))
    chr = strtoi(strs[1])
    site = strtoi(strs[2])
    locus = geno[,l]
    inc <- FALSE
    if (hetLocus(locus)) {
      num_het_loci <- num_het_loci + 1
      fit_pre <- fitnessFunc(meanG(pop))
      pop0 <- editGenome(pop, ind=c(1:popSize), chr=chr, segSites=site, allele=0, simParam=SP)
      pop1 <- editGenome(pop, ind=c(1:popSize), chr=chr, segSites=site, allele=1, simParam=SP)
      fit_post0 <- fitnessFunc(meanG(pop0))
      fit_post1 <- fitnessFunc(meanG(pop1))
      if ((fit_post0 > fit_pre) || (fit_post1 > fit_pre) ) {
        inc <- TRUE
      }
    }
    if (inc == TRUE) {
      num_inc <- num_inc + 1
    }
    if (FALSE) {
      for (i in 1:nInd(pop)) {
        ind = geno[i,]
        fit_pre <- fitnessFunc(pop@gv[i])
        pop0 <- editGenome(pop, ind=i, chr=chr, segSites=site, allele=0, simParam=SP)
        pop1 <- editGenome(pop, ind=i, chr=chr, segSites=site, allele=1, simParam=SP)
        
        fit_post0 <- fitnessFunc(pop0@gv[i])
        fit_post1 <- fitnessFunc(pop1@gv[i])
        if ((fit_post0 > fit_pre)) {
          inc <- TRUE
          break
        }
        if ((fit_post1 > fit_pre)) {
          inc <- TRUE
          break
        }
        
      }
    }
    
  }
  
  if (num_het_loci == 0) {
    perc <- 0
  } else {
    perc <- (num_inc / num_het_loci)
  }
  df$segAlleles[gen] <- num_het_loci
  df$percFitInc[gen] <- perc
  pop <- selectCross(pop=pop, trait=oneTraitFitFunc, nInd=popSize*(1-selIntensity), nCrosses=popSize)
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20

par(mar = c(5, 5, 3, 5) + 0.3)
plot(df$gen, df$fitness, type="l", lwd = "3", col=2, xlab = "Generation", ylab = "Fitness")
par(new = TRUE) 

plot(df$gen, df$percFitInc, type="l", lwd = "3", col = 3, axes = FALSE, xlab = "", ylab = "") 
axis(side = 4, at = pretty(range(df$percFitInc)))

mtext("% of Segregating alleles that are Beneficial", side = 4, line = 3)
par(xpd=TRUE)

legend("right",
  c("Trait Value", "% Alleles"),
       lty = 1,
       lwd = 3,
       col = 2:3)

ggplot(data=df, aes(x=gen, y=segAlleles)) +
  geom_line() +
  labs(x = "Generation", y = "Num Segregating Alleles")

Simulate several adaptive walks with different population sizes

nGens = 50
popSizes = c(10,200,1000)
nChr = 10
nSegSites = 50
selIntensity = 0.5
nLoci = nSegSites * nChr

fig <- plot_ly()
fit_df <- data.frame(gen=1:nGens,
                   fitness=numeric(nGens),
                   traitValA=numeric(nGens),
                   traitValB=numeric(nGens))
for (p in c(1:length(popSizes))) {
  popSize = popSizes[p]
  founders = runMacs(
     nInd=popSize,
     nChr=nChr,
     segSites=nSegSites
  )
  SP <- SimParam$new(founders)
  SP$addTraitA(mean=c(-10,10), var=c(1,1), nQtlPerChr=10)
  
  SP$setVarE(H2=c(0.5,0.5))
  foundingPop <- newPop(founders, simParam = SP)
  
  pop <- foundingPop
  
  for(gen in 1:nGens) {
    fit_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(pop)))
    fit_df$traitValA[gen] <- meanP(pop)[1]
    fit_df$traitValB[gen] <- meanP(pop)[2]
    pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=popSize*(1-selIntensity), nCrosses=popSize)
  }
  
  fig <- add_trace(
    fig,
    fit_df,
    name = popSize,
    x = fit_df$traitValA,
    y = fit_df$traitValB,
    z = fit_df$fitness,
    type = 'scatter3d',
    mode = 'lines',
    opacity = 1,
    color = p,
    line = list(width = 10)
  )
  
}

fig %>%
  layout(legend=list(title=list(text='Population Size')),
         scene = list(xaxis = list(title = "Trait A"),
                      yaxis = list(title = "Trait B"),
                      zaxis = list(title = "w"),
                      aspectmode='cube')) %>% hide_colorbar()

NA

Overlay an adaptive walk over a fitness landscape

set.seed(123)

nGens = 50
popSize = 200
nChr = 10
nSegSites = 50
nLoci = nSegSites * nChr

fig <- plot_ly()
fit_df <- data.frame(gen=1:nGens,
                   fitness=numeric(nGens),
                   traitValA=numeric(nGens),
                   traitValB=numeric(nGens))
founders = runMacs(
   nInd=popSize,
   nChr=nChr,
   segSites=nSegSites
)
SP <- SimParam$new(founders)
SP$addTraitA(mean=c(runif(1,-10,10),runif(1,-10,10)), var=c(1,1), nQtlPerChr=10)

SP$setVarE(H2=c(0.5,0.5))
foundingPop <- newPop(founders, simParam = SP)

pop <- foundingPop

for(gen in 1:nGens) {
  fit_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(pop)))
  fit_df$traitValA[gen] <- meanP(pop)[1]
  fit_df$traitValB[gen] <- meanP(pop)[2]
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=popSize*0.5, nCrosses=popSize)
}

overlayWalkOnLandscape(fit_df, type="SURFACE", calculateFitnessTwoTrait)

TODO:

Simualate mutations for a single individual to figure out P(allelic substitution is favorable) - should reflect Fisher/Kimura

P(mutation gets fixed) - show that this matches Kimura

Compare adaptive walks for DH/inbred population where there are new mutations VS population w/standing genetic variation

Figure out how to do QTL mapping

nBurnInGens <- 20

fit_df <- data.frame(gen=1:nGens,
                   fitness=numeric(nGens),
                   traitValA=numeric(nGens),
                   traitValB=numeric(nGens))

pop <- foundingPop

# BURN-IN
for (gen in 1:nBurnInGens) {
  fit_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(pop)))
  fit_df$traitValA[gen] <- meanP(pop)[1]
  fit_df$traitValB[gen] <- meanP(pop)[2]
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=popSize*0.5, nCrosses=popSize)
}

# Create 2 sub populations

popA_df <- data.frame(gen=1:nGens,
                   fitness=numeric(nGens),
                   traitValA=numeric(nGens),
                   traitValB=numeric(nGens))

popB_df <- data.frame(gen=1:nGens,
                   fitness=numeric(nGens),
                   traitValA=numeric(nGens),
                   traitValB=numeric(nGens))

# Figure out a better way to select randomly? or by traitA, traitB?
popA <- selectInd(pop, use="rand", nInd=popSize*0.1)
popB <- selectInd(pop, use="rand", nInd=popSize*0.1)

nSegGens <- 20
for (gen in 1:nSegGens) {
  popA_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(popA)))
  popA_df$traitValA[gen] <- meanP(popA)[1]
  popA_df$traitValB[gen] <- meanP(popA)[2]
  popA <- selectCross(popA, trait=twoTraitFitFunc, nInd=nInd(popA)*0.5, nCrosses=nInd(popA))
  
  popB_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(popB)))
  popB_df$traitValA[gen] <- meanP(popB)[1]
  popB_df$traitValB[gen] <- meanP(popB)[2]
  popB <- selectCross(popB, trait=twoTraitFitFunc, nInd=nInd(popB)*0.5, nCrosses=nInd(popB))
}

Pull out individuals from each population and cross them Create RILs Inspect genetic architecture + do QTL mapping Work in Progress

Figure out QTL Mapping